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ABSTRACT 

Context. Numerical simulations of stellar convection and photospheres have been developed to the point where detailed shapes of 
observed spectral lines can be explained. Stellar atmospheres are very complex, and very different physical regimes are present in the 
convection zone, photosphere, chromosphere, transition region and corona. To understand the details of the atmosphere it is necessary 
to simulate the whole atmosphere since the different layers interact strongly. These physical regimes are very diverse and it takes a 
highly efficient massively parallel numerical code to solve the associated equations. 

Aims. The design, implementation and validation of the massively parallel numerical code Bifrost for simulating stellar atmospheres 
from the convection zone to the corona. 

Methods. The code is subjected to a number of validation tests, among them the Sod shock tube test, the Orzag-Tang colliding shock 
test, boundary condition tests and tests of how the code treats magnetic field advection, chromospheric radiation, radiative transfer in 
an isothermal scattering atmosphere, hydrogen ionization and thermal conduction. 

Results. Bifrost completes the tests with good results and shows near linear efficiency scaling to thousands of computing cores. 
Key words. Magnetohydrodynamics (MHD) - Radiative transfer - Methods: numerical - Sun: atmosphere - Stars: atmospheres 



1. Introduction 



and heating of the chromosphere and corona (Abbett 2007 



The development of faster computers and better algorithms has 
made simulations a viable experimental tool to understand a 
number of astrophysical problems in detail. This development 
was clear more than a decade ago (Miyama et al. 1999 1. The 
years that have followed have borne this out fully and there are 
now a number of groups that are modeling the outer layers of 
cool stars, including magnetic fields, to a high degree of realism 



CGudiksen & Nordlund'2005','Hansteen'20 04t|Steiii & Nordlund 
2006, Hansteen et al. 2007; Schaffenber ger et al.||200 5; Vogler 
et al.l|2005t |Heinemann et al.i|2007i |Abbett||2007nTsobe et al. 
20081 1. 

Vital to this effort was the development of multi-group tech- 
niques to handle radiative transfer in the photosphere (Nordlund] 
1 1982) and for mode ls extending into the chromosphere, scatter- 
ing ( Skartlien 2000 1. The majority of the codes mentioned above 
employ these multi-group methods. 

As the initial exploratory phase for codes including mag- 
netoconvection is nearing successful completion, a number of 
challenging problems are now being considered with some con- 
fidence. 

The problems include issues such as the existence and for- 
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Recent codes solve the full radiative magnetohydrodynamic 
(MHD) equations with some precision. However, additional 
physics may have to be considered in order to solve problems 
inherent to the low density, low and/or high temperature condi- 
tions of the outer solar atmosphere. These effects encompass the 



inclusion of thermal conduction by various methods ( [Gudiksen 
|& Nordlund||2005[ [Hansteen et al.||2007| l, but also GeneraUzed 



Ohm's law ( Arber et al. 2007[ l, aswell as non-equilibrium effects 



|2009a[ [Tortosa-Andreu & Moreno-InsertiS||2009^ the structure 



such as hydrogen ionization and molecule formation ( Leenaarts 
|et al. 2007; Wedemeyer-Bohm & Steffen 2007 ), are now being 
added to these codes. 

Numerical simulations require complex algorithms to solve 
the physics required, but in addition, combining high spatial res- 
olution with the large spatial scales characteristic of atmospheric 
phenomena requires large memory and many CPU hours com- 
puting time. The cost of high performance computing (HPC) 
specific hardware has driven the market for supercomputers to be 
focused mainly on utilizing relatively cheap off-the-shelf com- 
puter parts instead of developing specific supercomputing hard- 
ware. The cheapest performance can be gained from connecting 
almost standard computers through a local network. The inter- 
connect speed and communication software implementation is 
what sets one cluster apart from another. Such clusters have a 
distributed memory architecture, which is different from previ- 
ous generations of supercomputers with specifically built hard- 
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ware which were typically shared memory vector machines. The 
consequence is that the numerical codes used in for instance as- 
trophysics have to be highly parallelized in order to perform well 
on large HPC systems using, for example, the Message Passing 
Interface (MPI). In the near future we expect further develop- 
ments of these and other techniques such as the widespread 
use of adaptive mesh refinement and graphics processor units 
(GPUs). 

In this paper we present a code, BifrosiQ to solve the MHD 
equations in a stellar atmosphere context, specifically designed 
to take advantage of the environment provided by modern mas- 
sively parallel machines connected through potentially slow 
communication channels. The main design goal has been to re- 
duce the amount of communication required at each timestep 
while retaining good scaling on problems requiring up to several 
thousand processors. 

In addition, we aim to create a flexible design in which var- 
ious physical processes and extensions to the basic MHD equa- 
tions are straightforward to implement without rewriting large 
portions of the code. 



2. Basic equations and set-up 

Bifrost is built on several generations of previous numerical 
codes, where the Oslo Stagger Code is the latest, which at 
their core have the same implemented method (Nordlund & 
|Galsgaard|1995| palsga ard & Nordlun d 1996). The core of the 
code remains the same, but as the previous generations of codes 
required shared memory architectures, the need for a new mas- 
sively parallel code able to run on distributed memory comput- 
ers was obvious. As the core of the code was rewritten, we had 
the opportunity to make Bifrost much more modular and user 
friendly than the previous generations of codes. 

At the core of the code is a staggered mesh explicit code that 
solves the standard MHD partial differential equations (PDEs) 
on a Cartesian grid: 



dpu 

-h- = -V ■ (puu - T)-VP + J xB + pg 
at 

E = TjJ -uxB 

dB „ ^ 
— = -VxE 

dt 
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dt ^ 



(1) 

(2) 

(3) 
(4) 

(5) 
(6) 



where p,u,e,B are the density, the velocity vector, the internal 
energy per unit volume and the magnetic flux density vector re- 
spectively. T ,P,J, g,fi,E and are the stress tensor, the gas 
pressure, the electric current density vector, the gravitational ac- 
celeration, the vacuum permeability, the electric field vector and 
the magnetic diffusivity respectively. The quantity Q can con- 
tain a number of terms, depending on the individual experiment. 
It could for instance contain a term from the chosen Equation Of 
State(EOS), a term containing the effect of the Spitzer thermal 
conductivity, a term from radiative transfer, etc. The EOS needed 
to close this set of equation can be anything from a simple ideal 



' In Norse mythology, Bifrost (pron. bee-frost) is the name of the 
rainbow bridge from Midgard (the real of man) to Asgard (the realm of 
the gods), build from fire, water and air. 



gas EOS to a complex EOS including detailed microphysics (See 
Sect. IT). 

3. The method 

The basic operator in the code is the 6th order differential oper- 
ator. The derivative of the function / in the grid point location 
+ 1/2 has the form 

dx 
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where the subscript is the gridpoint number and a^, b^, depend 
on the grid point distance. As the derivative operator produces 
results that lie half way between grid points, the variables are not 
co-located in space, but placed on staggered grids such that some 
variables are placed at cell centers, some on cell faces and some 
at cell corners. By carefully choosing which variables to place 
in each of these locations, it is possible to minimize the number 
of interpolations needed in order to realign computed values in 
space. Nevertheless, it is not possible to escape interpolations, 
and when it is needed the interpolation procedures used are 5th 
order and look very similar to the derivative operators: 

/+1/2 = a[/o + 
b[f-i+U2] + 

c[/-2+A3] (8) 

where a, b and c are constants. 

The computational grid can be stretched in one direction at a 
time, meaning that dx is not constant in the computational vol- 
ume. The stretching of the grid is done by a simple Jacobian 
transformation after a derivative operator is used. Stretching in 
more than one direction cannot be accomplished in this simple 
manner without using a number of interpolations which would 
decrease precision, and at the moment it is not possible to em- 
ploy an adaptive mesh refinement scheme in the code, but due 
to the modularity of the code this can be implemented at a later 
stage. 

3.1. Diffusion 

All numerical codes are diffusive in nature, just due to the dis- 
crete nature of the algorithms used in solving the equations and 
because of the inaccurate machines used to solve them on. This 
is even true for implicit codes that do not contain any direct dif- 
fusive terms in their equations, that said, the diffusion in implicit 
codes is significantly smaller than for explicitly diffusive codes 
of the same order. Bifrost is an explicit code and it is therefore 
necessary to include diffusive terms in the Eq. [TJ|6] in order to 
maintain stability. Bifrost employs a diffusive operator which 
is split in two major parts: A small global diffusive term and a 
location-specific diffusion term (sometimes dubbed "hyper dif- 
fusion"). The diffusive operator in one dimension is of the form 



dt ■ 
where 
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dx 
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and Vj is the first order gradient in the x direction. The splitting 
of the diffusive terms into local and global components makes it 
possible to run the code with a global diffusivity that is at least 
a factor 10 less than if the global term were the only one im- 
plemented in the code. The splitting also makes it impossible to 
provide a single Reynolds number, Magnetic Reynolds number 
or any other dimensionless number that includes the diffusion 
constant, since diffusion is not constant in space or time when 
running Bifrost. Therefore it is only possible to provide a range 
of the above mentioned numbers, of which the smaller, global, 
value for the diffusion would be correct for most of the simu- 
lation volume most of the time (unless, for instance, simulating 
supersonic turbulence), while the higher number of the range for 
the diffusion would be valid in locations characterized by very 
large values of the gradients of a certain variable. Thus, in prin- 
ciple, the code can be run at much greater values of the relevant 
Reynolds numbers in most of the computational domain than 
would be otherwise possible using a single diffusion coefficient. 



4. Modules 

Bifrost has been created with a high degree of modularity. The 
code has a basic skeleton which connects to a number of mod- 
ules. Any of these modules can contain a simple procedure or 
method, include a number of sub-modules or be left empty. For 
example, the timestepping procedure can be swapped between a 
Runge-Kutta scheme and a Hyman scheme by changing only a 
single line in an input file to the compiler, plus a recompilation. 
Most interestingly, the code can include any number of modules 
that provide new physics, or boundary conditions, in the same 
simple way. 

As there can be several implemented modules handling the 
same job in the code (timestepping, EOS, radiative transfer etc) 
this section contains the most important of these, which are pre- 
sented either with a short description if the modules are stan- 
dard numerical schemes, or in more detail if they are specific to 
Bifrost. 



5. Timestepping 

Timestepping can be handled by two different procedures: 
A third-order Runge-Kutta method or a third-order Hyman 
method. Both of these procedures produce nearly identical re- 
sults. The Runge-Kutta method is able to handle a longer 
timestep, while being more computationally intensive, so in 
terms of CPU time their effectiveness is almost the same. 



5. 1 . Third-order Runge-Kutta timestepping sclieme 

The implemented 3rd order Runge-Kutta scheme splits the 
timestep into three sub-steps. In order to take one timestep all 
the partial differential equations are solved three times but it 
is not necessary to save more than two results at a time in the 
same timestep, making this scheme 2N in memory requirements 
(where in the amount of memory needed to run through the 
partial differential equations once). For large simulations that 
can be a considerable memory saving feature. The advantage 
in this method is that the intermediate timestep results can be 
used to extrapolate the result of the total timestep further in time 
than that possible by using three separate time steps, while at 
the same time having a high order precision. The Runge-Kutta 



method is defined by assuming the change in the variable / can 
be written 



(11) 



The change during one timestep At of the variable / is then given 
by 



Af 

fit + At) = fit) + (^i + 4k2 + h) , 
o 

where 

ki = Fit, fit)) , 

/I 1 \ 

k2 - Flt+-At,fit)-\--kiAt\ and 

hi = Fit + AtJit)-kiAt + 2k2At) . 



(12) 

(13) 
(14) 
(15) 



5.2. Third-order Hyman timestepping scheme 

The third-order Hyman predictor-corrector scheme is described 



by Hyman ( 1979 1. It is an iterative multistep method employ- 



ing a leap-frog scheme to attain 3rd order accuracy in time. The 
method is quite simple and can be described by assuming that 
the differential equations can be written in the following form: 



Fit Jit)) 



(16) 



The Hyman method's first step is to find a second order predic- 
tive solution to Eq. [T6]by using the formula: 



r. 



<2) 



il - r^)f„ + rj„^i + Atil + r)F„ 



(17) 



where the superscript is the order of the term, the subscript is the 
timestep number and F it, /„(?)) has been written F„ for simplic- 
ity and 



r = (Af„+i/Af„) 



(18) 



with Af being the timestep. Then the PDE's are solved again and 
finally a corrector is applied given by 



= [i2-r)il+r)^f„ + r'fn-i + 

Af(l + r)V„ + Af(l + r)Ff^\] /(2 + 3r) 

6. Boundary conditions 



(19) 



Bifrost is able to make the computational box periodic in one, 
two or three dimensions. Non-periodic boundaries imply that a 
boundary condition must be imposed on that boundary. Bifrost 
implements non-periodic boundary conditions by padding the 
computational domain with ghost zones in the relevant direction 
and filling them according to the boundary condition chosen. 
Several standard boundary conditions are implemented includ- 
ing symmetric, antisymmetric as well as extrapolated boundary 
values that fill the ghost zones before the MHD PDEs are com- 
puted. Experiment-specific boundary conditions such as con- 
stant entropy flux for stellar convection simulations and char- 
acteristic boundary conditions can also be used by adhering to a 
simple format for boundary calls. In the latter case of character- 
istic boundaries the conservative MHD PDEs are replaced by the 
characteristic PDEs at the boundary zones before the variables in 
these zones are updated in the usual manner. 
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6.1. Characteristic boundary conditions 



7. Equation of state 



The characteristic boundary conditions aim at transferring dis- 
turbances through the boundaries without any or at least with 
minimal reflection, a problem that can plague standard bound- 
ary conditions that simply use symmetric or extrapolated values. 
In this case the 8 MHD equations are written in terms of the 
characteristics and split into horizontal and vertical components 
in the following manner. 

With U defining a vector that contains the conserved MHD 
variables we may write the equations as 



5U' ^^ ^, 
+ / -7. — = D , 

at ^ Bx,. 



k=\ 



(20) 



where F contains the fluxes, D the source terms, and m - 3 for 
a 3D problem. These conservation equations can be transformed 
using linear algebra into "primitive", wave-like equations for a 
corresponding set of 8 primitive variables U, 



at £^ dxk 



(21) 



where the three A'^ are 8 x 8 matrices. The choice of primitive 
variables U is not unique, and in principle could be taken to be 
the conserved variables. The goal of this procedure is then to ar- 
rive at equations that resemble simple advection equations where 
specifying boundary conditions is straightforward: extrapolation 
based on one-sided derivatives for outgoing characteristics and 
on the basis of exterior data, such as no incident waves, for 
incoming characteristics. By combining all the flux divergence 
terms except those in the direction perpendicular to the bound- 
ary, now called z, with the source term (forming a new source 
term C) we can write characteristic equations for the z direction 
in the sought after form 



ot oz 



(22) 



as shown by e.g. Thompson ( 1987| l. Here, the rows ; of the ma- 
trix S"' are given by the left eigenvectors 1^, and A is the diago- 
nal matrix formed by the eigenvalues /I, of the matrix A*^ belong- 
ing to the z dkection. We now define a vector d containing the z 
derivatives of the characteristic equations as 

OZ 

Having isolated the various characteristic wave modes propagat- 
ing in the z direction we left-multiply Eq. 22 by S in order to 



write the MHD equations in primitive form in terms of d: 



dV 

~dt 



-HSd^C 



(23) 



where the primitive variables U are comprised of the variables 
p, u, e and B. The d vector containing the z derivatives di - d% 
constitute the information that is flowing along the characteris- 
tics; outflowing characteristics are defined by the interior solu- 
tion, inflowing characteristics by the requirement that the incom- 
ing wave be constant in time. Expressions for Eq. 23 in primi- 



tive form and for the characteristic derivatives d can be found in 
Appendix [A| 



The code provides several different EOS modules, which can be 
chosen according to the experiment one wants to perform. The 
EOS modules provide the temperature and pressure and their 
thermodynamic derivatives as a function of mass density and in- 
ternal energy per mass. There are currently three different mod- 
ules. 

The first implements an ideal gas EOS, suited for testing and 
idealized experiments. 

The second module implements an EOS based on tables 
generated with the Uppsala Opacity Package (Gustafsson et al. 



19751. It assumes local thermodynamic equilibrium (LTE) for 



atomic level populations and instantaneous molecular dissocia- 
tion equilibria. This package is required when running with full 
radiative transfer (see Sect. 8.3 i, as it also provides the opac- 
ity, thermal emission and scattering probability for the radiation 
bins. The tables are generated with a separate program; different 
tables can be generated to account for different stellar spectral 
type, chemical composition and number of radiation bins. 

The third package computes the gas temperature, gas pres- 
sure and electron density explicitly based on the non-equilibrium 
ionization of hydrogen in the solar atmosphere. This package 
can only be used for simulations of the solar atmosphere, as it 
depends on a number of parameters that vary with stellar spec- 
tral type. These parameters have so far only been determined for 
the sun. More details on non-equilibrium hydrogen ionization 
are given in Sect. |9] 

8. Radiative transfer 

Full 3D radiative transfer is computationally very costly. The 
properties of radiation are very different from the strictly local 
problem of MHD, since radiation can couple thermodynamic 
properties over arbitrarily long distances. For a code that relies 
on the local properties of MHD to parallelize, this makes radia- 
tive transfer costly not just in shear computations needed to solve 
the problem, but also because the results of the computation must 
be communicated to all nodes. As a result, the modules handling 
radiative transfer in Bifrost have employed assumptions which 
simplify the calculations to some extent. At the moment three 
modules handle radiative transfer depending on the problem at 
hand and they are often combined. 

The electron number density is often needed in the radia- 
tive transfer modules (e.g. for the calculation of collisional ex- 
citation rates and opacities). The electron number density com- 
putation depends on which EOS package is used: if hydrogen 
non-equilibrium ionization is calculated, the electron number 
density comes from the simultaneous solution of the hydrogen 
rate equations, the energy equation and the charge conservation 
equation (see Sect.|9]for details). For the EOS package based on 
tables generated with the Uppsala Opacity Package, the electron 
density comes from solving the Saha equation for all species. If 
the electron density is not provided by the EOS package but is 
needed (e.g., for the chromospheric radiation approximation, see 
Sect. |8.2| l, it is computed using the Saha relation for hydrogen, 
but setting a floor to the ionization degree of 10""^ to account for 
the easily ionized metals in the solar atmosphere. 



8.1. Optically thin radiative transfer 

In the outer solar atmosphere (and many outer stellar atmo- 
spheres) radiative losses can be simply treated assuming that 
the atmosphere is optically thin. In the sun this is true for most 
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lines from the upper chromosphere/lower transition region up 
to the corona. In this case (and assuming ionization equilibrium 
only dependent on temperature) the radiative transfer problem 
reduces to a radiative loss which only depends on density and 
temperature of the form (note that positive Q corresponds to 
heating): 



Gthin = -ne nJ{T) exp(-P/Po) 



(24) 



where riu and Wg are the number densities of hydrogen and elec- 
trons respectively, f(T) is a function of the temperature and the 
term exp(-P/Po) provides a cutoff where P > Pq. The radiation 
stems mainly from the resonance lines of multiply ionized ele- 
ments such as carbon, oxygen, and iron and the function /(T) 
can be pre-computed assuming ionization equilibrium. This re- 
lation is valid from roughly 2 x 10'* K and up. In Bifrost the 
function f{T) is computed by using ionization and recombina- 
tion rates given by |Arnaud & R othenflug (1985j) and jShuU &| 
van SteenbergI ( |1982| l and collisional excitation rates given by 
the HAO-DIAPER atom data package ( [Judge & Meisnerj 1994) 1 
including lines from He, C, O, Ne and Fe. The optically thin 
losses from hydrogen lines are normally calculated in the chro- 
mospheric approximate radiation part, see Sect. 8.2 but may in- 
stead be included here in the cases where a calibration of the 
chromospheric radiative losses is lacking. The total hydrogen 
number density is derived from the plasma density p assum- 
ing solar abundances. We set Pq to a typical value of the mid- 
chromospheric pressure and the term exp(-P/Po) then makes 
sure that there are no optically thin contributions calculated in 
the deep convection zone where the temperature is also above 
2 X lO-* K. 



8.2. Chromospheric radiation approximation 

Chromospheric radiative losses are dominated by a small num- 
ber of strong lines from hydrogen, calcium and magnesium 
The source functions and opacities of 



(Vemazza et al. 1981 



these lines are very much out of local equilibrium and the op- 
tical depth is significant; neither the optically thin approach of 



Sect. 8.1 nor the full radiative transfer with coherent scattering 
and LTE opacities (see Sect. 8.3 i give good results. We approxi- 
mate the radiative loss in these lines with the formulae 



2[H,Ca,Mg] - -C(r)[H,Ca,Mg]neP<^[H,Ca,Mg]('«c) 



(25) 



where C(T)nsp gives the total collisional excitation rate (in en- 
ergy per volume per unit time), (pim^.) gives the probability that 
this energy escapes from the atmosphere and is the column 
mass. The function 0(mc) and the temperature-dependent co- 
efficient C{T) are determined for each element from detailed 
ID radiative transfer computations with the RADYN code (see, 
e.g.,|Carlsson & Stein|l9 95 ) and 2D computations with Multi3D 
( Leenaarts & Carlsson|2009) . These functions include hydrogen 
lines and the Lyman continuum and all lines and continua from 
Ca II and Mg II. The method is described in detail in Carlsson & 
Leenaarts (in preparation). 

Half of the UV radiation lost from the corona in optically 
thin lines (see Sect. |8.1| l goes towards the sun and most of that is 
eventually absorbed in the chromosphere, mostly in the Lyman 
continuum and the He I continuum. This radiative heating is 
modeled through the representative bound-free absorption cross- 
section cr of He I at 25 nm: 



Qhs - O" He,25nm "He I eXp(-THe,25 nm) Jt 



thin 



with /thin the angle-averaged radiation field caused by Qthin (see 
Carlsson & Leenaarts in preparation). 



8.3. Full radiative transfer 

Full radiative transfer computations are required when a simula- 
tion includes the convection zone beneath the photosphere, cov- 
ering optically thick regions, optically thin regions, and the tran- 
sition between the two regimes. A simplified treatment using, 
e.g., Newtonian cooling or the diffusion approximation, cannot 
provide sufficiently realistic radiative heating and cooling rates 
in this boundary layer. 

Owing to the very short time scales of photon interaction and 
propagation in a convective stellar atmosphere, it is possible to 
solve radiative transfer as a time-independent problem, resulting 
in the expression 



n ■ V/^(x, n) = -;t'.i(x)/^(x, n) -i- j^ix, n) . 



(27) 



where /^(x, n) is the monochromatic specific intensity at spatial 
point X for a beam in direction n at wavelength A, xi is the gas 
opacity, and j,\ is the emissivity. The two material constants xa 
and depend both on the thermodynamic state of the gas and 
on the radiation field, and are highly wavelength-dependent in 
the presence of spectral lines and other atomic and molecular 
transitions. Velocity fields, such as convective motions in the at- 
mosphere, lead to an additional coupling between ray directions 
and wavelengths through Doppler shifts. The complexity of tak- 
ing full account of the underlying physical mechanisms is com- 
putationally prohibitive. We therefore assume a static medium 
and LTE for the gas opacity, which thus depends only on the 
local gas temperature and the gas density, and which can there- 
fore be precomputed and tabulated. [Skarthen ( 2000[ l and Hayek 
et al. (2010j) showed the importance of photon scattering in sim- 
ulations of the higher solar atmosphere; we therefore include a 
coherent scattering contribution in the gas emission, which re- 
quires an iterative solution of the radiative transfer equation to 
obtain a consistent radiation field. 

The vast number of spectral lines encountered in a stellar 
atmosphere requires, in principle, the solution of a large num- 
ber of radiative transfer problems. This is currently too de- 
manding for realistic 3D radiation-hydrodynamical calculations. 
We therefore approximate the opacity spectrum by substitut- 
ing the monochromatic xa with a small number of mean opac- 
ities (Nord lund 1982; Skartlien 2000) and solving wavelength- 
integrated radiative transfer. 

The radiation field encountered in cool stellar atmospheres 
does not contribute significantly to the momentum balance. We 
consider only a radiative heating rate Qiadj for every mean opac- 
ity (index /), given by the first moment of Eq. [27] 



a-ad,; = -V ■¥i^'Xnxi{Ji-S[) . 



(28) 



where Fi is the radiative energy flux, 7, is the mean intensity 
and Si = jilxi is the source function. The solver computes the 
radiation field using the short characteristics technique in every 
spatial subdomain and iterates the solution until convergence. 
The radiative heating rates are then added to the energy equation 
(Eq. |6]l. If the hydrodynamical timesteps are sufliciently short 
that changes in the radiative energy flux are small, full radiative 
transfer only needs to be computed in a fraction of the hydrody- 
namical timesteps. A detailed description of the numerical meth- 
ods and the parallelization techniques used for the full radiative 
transfer module in Bifrost can be found in [Hayek et al.| ( )2010f 



(26) 9. Hydrogen ionization 



The ionization of hydrogen in the solar chromosphere does not 
obey LTE or statistical equilibrium ( [Carlsson & Stein]|2002| l. 
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Proper modeling of this ionization requires that non-equilibrium 
effects be taken into account. 

The hydrogen ionization module solves the time-dependent 
rate equations for the atomic hydrogen level populations «, 



dn 



(29) 



together with an equation for time-dependent H2 molecule for- 
mation and equations for energy and charge conservation. Here 
Pij is the rate coefficient for the transition from level / to j and n\ 
is the number of energy levels in the hydrogen atom (normally 
set to 6 in Bifrost). Solution of this system of equations yields 
the gas temperature, electron density and the atomic and molec- 
ular hydrogen populations (nH2)- The radiative transition rate co- 
efficients in the rate equations are prescribed as determined by 
|Sollum| ( |1999| l. This removes the effect of the global coupling of 
the radiation and makes the problem computationally tractable, 
but still computationally demanding. 
The gas pressure is computed as 



P = kT 



"I 



(30) 



with «o the number density of all other atoms and molecules that 
are not, or do not contain, hydrogen. 

Full non-LTE radiation tables as functions of and T can be 
used to replace the radiation tables as functions of e and p, which 
assume hydrogen ionization based on Saha ionization equilib- 
rium. 

The method is described in detail in |Leenaarts et al.| ( |2007[ l, 
the additional equation for time-dependent H2 formation is given 
in Appendix [B] 

10. Thermal conduction 

As the plasma temperature rises towards one million degrees in 
the tenuous transition region and corona, thermal conduction be- 
comes one of the major terms in the energy equation, and mod- 
eling of this term is vital if this portion of the atmosphere is to be 
simulated with any fidelity. The implemented form of the ther- 
mal conduction takes the form ( |Spitzer|1956| ) : 



(31) 



where the gradient of T is taken only along the magnetic field 
(V||) and kq is the thermal conduction coefficient. The conduc- 
tion across the field is significantly smaller under the conditions 
present in the solar atmosphere and is smaller than the numerical 
diffusion, so it is ignored. Since thermal conduction is described 
by a second order diffusion operator, this introduces several dif- 
ficulties: The Courant condition for a diffusive operator such as 
that scales with the grid size Ax^ instead of with Ax for the mag- 
netohydrodynamic operators. This severely limits the time step 
Af the code can be stably run at. We have implemented two so- 
lutions to this problem. 

In the first, the thermal flux is calculated and if the diver- 
gence of the thermal flux sets severe restraints on the timestep, 
it is throttled back by locally lowering the thermal conduction 
coefficient kq. This method is only acceptable when the number 
of points where the conduction has to be throttled back is very 
small and the results must be analyzed carefully. 

A second solution is to proceed by operator splitting, such 
that the operator advancing the variables in time is L = ihydro + 



^conduction- The conductivc part of the energy equation is handled 
by discretizing 



5e 
di 



-v-F, = -V||-[^or'^'V||r] 



(32) 



and solving the resulting problem implicitly. 

We discretize the Lc„„d,iciion operator using the Crank- 
Nicholson (or '0') method and thus write 



(e«+' - e*)/dt = 0V|| ■ + (1 - 0)V|| ■ F«) 



(33) 



where the quantities with superscript n are computed before the 
hydrodynamic timestep, the quantities with superscript * are 
computed after the hydrodynamic timestep and the quantities 
with superscript n + I are the temperatures deduced implicitly. 
The variable is set to a value between 0.5 and 1 . The implicit 
part of the problem is in our case computed using a multi-grid 
solver (A concise introduction to multi-grid methods may be 
found in chapter 19 of Press et aL] ( |1992[ )). 

In implementing the multi-grid solver, we use the same do- 
main decomposition as in the magnetohydrodynamic part of the 
code. The small scale residuals in Eq. 33 are smoothed using a 
few Gauss-Seidel sweeps at high resolution, before the result- 
ing temporary solution is injected onto a coarser grid on which 
smoothing is again performed. The process continues by Gauss- 
Seidel smoothing sweeps on steadily coarser grids until the size 
of the problem on individual processors is small enough to be 
communicated and stored on each and every processor. At this 
point the partial solutions are spread to all processors which con- 
tinue the Gauss-Seidel smoothing, now on the global problem, 
on steadily coarser grids. Finally, when the coarsest grid size is 
reached, the solution is driven to convergence by performing a 
great number of Gauss-Seidel sweeps. Subsequently, the coarser 
solutions are corrected and then prolonged on successively finer 
grids, first globally, but after a certain grid size is reached the 
temporary solution is spread and the problem is again solved lo- 
cally, using the same domain decomposition as in the magneto- 
hydrodynamic part of the code. After each prolongation, the so- 
lution is smoothed using a few Gauss-Seidel sweeps. The entire 
cycle can be repeated as many times as desired, until a converged 
solution is found. 



1 1 . Tests and validation 

Bifrost has been extensively tested in order to validate the results 
it provides. When possible, the tests are taken from the literature 
to validate the code and allow comparison with previous work. 
These tests were selected for their simplicity, their utility and 
their challenge to the different terms in the algorithms. 

11.1. The Sod shock tube test 



The Sod shock tube test ( |Sod|1978] l has become a standard test 
in computational HD and MHD. It consists of a one-dimensional 
flow discontinuity problem which provides a good test of a com- 
pressible code's ability to capture shocks and contact discontinu- 
ities within a small number of zones, and to produce the correct 
density profile in a rarefaction wave. The test can also be used to 
check if the code is able to satisfy the Rankine-Hugoniot shock 
jump conditions, since this test has an analytical solution. 

The code is set to run a ID problem and using the initial 
conditions for the Sod problem. The fluid is initially at rest on 
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either side of a density and pressure jump. The density and pres- 
sure jumps are chosen so that all three types of flow discontinu- 
ity (shock, contact, and rarefaction) develop. To the left, respec- 
tively right side of the interface we have: 



PR = 1 , 

PL = 0.125 , 

Pr = 1/7, 

Pl = 0.125/r. 



(34) 
(35) 
(36) 
(37) 



The ratio of specific heats is chosen to be y = 5/3 on both 
sides of the interface and with a uniform magnetic field perpen- 
dicular to the ID axis of the domain {B, = 1). The units are 
normalized, with the density and pressure in units of the density 
and pressure on the left hand side of jump and the velocity in 
units of the sound speed. The length is in unit of the size of the 
domain and the time in units of the time required to cross the 
domain at the speed of sound. 

We have run the simulation at dififerent resolutions and with 
different values of the numerical diffusion coefficient in an effort 
to find the minimal diffusion the code can be run with without 
developing numerical instabilities. From this test and the follow- 
ing "field advection" and Orzag-Tang tests we find that the dif- 
fusion parameters should be greater than vi > 0.02, V2 > 0.2 and 
V3 > 0.2 to ensure stability. Other parameters, such as the quench 
parameter, have also been varied in testing the code, in order to 
find the best values for capturing the Sod shocks properly. 

Fig.[T]shows the density and pressure profiles at time 0.193. 
The solid line is our numerical solution and the dashed line is the 
analytic solution at the same instant in time. It is clear that when 
running with the numerical diffusion coefficients large enough 
to avoid post shock numerical instabilities, but low enough not 
to lose the sharp shock profiles, the code solves the Sod shock 
correctly. 

11.2. The magnetic field advection test 

This is a multidimensional convection test, which serves to test 
the conservation of the various MHD quantities such as the den- 
sity, momentum, and magnetic field with advection. Moreover, 
multidimensional MHD problems present a special challenge to 
the conservation of the divergence of the magnetic field (Toth 
|& OdstrciI||1996| l. This can only occur if the scheme preserves 
V • B = 0. 

This advection test is based on a test described previously by 
[DeVore (1991). The test involves advecting a cylindrical current 
distribution, which forms a tube of magnetic field, diagonally 
across the grid. Any arbitrary advection angle can be chosen and 
the tube can have any orientation. For the 3D results shown here, 
the problem domain is given by < .xi < 1; -1/(2 cos (30°)) < 
X2 < 1/(2 COS (30°)), < X3 < 1, and the flow is inclined at 
60 degrees to the X2 axis in the same plane as xi; where xi, X2 
and xj, could be either x, y or z. The loop is oriented in the x^ 
direction. This geometry ensures the flow does not cross the grid 
along a diagonal, so the fluxes in xi, X2 and X3-directions will be 
different. 

The flow velocity is set to 1.0, with direction ui = sin (60°) 
and M2 - cos (60°). We have run many tests, changing Ux, Uy and 
Mj between these ui and U2 values and have checked that there 
is no dependence on any direction. The density is 1, pressure is 
l/y, and the gas constant is y = 5/3. Periodic boundary condi- 
tions are used everywhere. 
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Fig. 1. Density profile (top panel) and pressure divided by the 
maximum of the pressure profile (bottom panel) at time 0.193. 
The solid line is the numerical result and the dashed line is the 
analytical solution. 



The magnetic field is initialized using a vector potential, to 
make sure that the magnetic field is initially divergence free. The 
potential in this example is given by: 



A3 = 0.009e"<^'' 

where r is the distance from the box centre in the x - y plane, 
and the two other components of the vector potential are set to 
zero. The magnetic vector potential provides the magnetic field 
vector according to B = V x A, and leaves a narrow cylinder of 
magnetic field and the divergence of the magnetic field is zero 
as calculated by Bifrost. 

The code has solved this test with different setups of the di- 
rection of the initial velocity and orientation of the loop. We have 
also checked that the test performs successfully with varying dif- 
fusion and quench parameters. The increase in V B with time lies 
within the numerical error. This means that every 10^ time-step 
the cumulative error of the numerical errors makes it necessary 
to perform a cleaning of the V ■ B. 
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Fig. 2. The magnitude of the magnetic field shown at f = (top) 
and at f = 1.155 (bottom) when the cylinder has moved one box 
length along the x-axis, and roughly half way along the y-axis. 
The lower contour plot has been centered to make a comparison 



In Fig. |2j the magnetic cylinder shows small alterations in 
shape due mainly to the rectangular grid used, and because the 
width of the magnetic ring initially was chosen to be at the limit 
of the resolution capability of the code, which exaggerates the 
effect of the grid. 

1 1.3. The Orszag-Tang test 



The problem was first studied by |Orszag & Tang] ( |1979| l. Since 
then solutions to the problem have been extensively compared 
in numerical MHD simulations. This provides a good way to 
compare codes. 

The Orszag-Tang (O-T) vertex is a model problem for test- 
ing the transition to supersonic 2D MHD turbulence. Thus, the 



problem tests how robust the code is at handling the forma- 
tion of MHD shocks and shock-shock interactions. The problem 
can also provide some quantitative estimates of how significant 
magnetic monopoles affect the numerical solutions, testing the 
V ■ B = condition. Finally, as mentioned, the problem is a very 
common test of numerical MHD codes in two dimensions, and 
has been used in many previous studies. 

The set up is the following: The domain is 2D and goes 
from < X < 1 and < y < 1. The boundary conditions 
are periodic. The density p = 1, the pressure P - 1/y and 
y = 5/3 everywhere. Note that this choice gives a sound speed 
of Cs = -iJyP/p - 1 • The initial velocities are periodic with: 



Ux — — sinilny) , 



(38) 
(39) 



The magnetic field is initialized using a periodic vector potential 
defined at zone corners: 



cos(4;7rx) cos(27r3;) 



4-n 



2n 



(40) 



with Bo - Vl/(4^)- Face-centered magnetic fields are computed 
using B = V X to guarantee V ■ B = initially. This gives: 



Br 



-Bo sin(27ry) . 
Bo sin(47rx) 



(41) 
(42) 



We have run a number of simulations with different grid 
resolution, diffusion and quench parameters. We observed good 
evolution of the shocks in all runs when we used diffusion pa- 
rameters set to vi > 0.02, V2 > 0.2 and V3 > 0.2 in agreement 
with the previous tests. However, the shocks will be smoother 
when increasing the diffusion or when decreasing the resolution. 
Fig. [3 is created to be directly compared with Fig. 3 of Ryu et al. 
P998 1 who use a upwind, total variation diminishing code, and 
even Ihough there are small differences it is hard to tell which 
code is superior in spite of the two codes being fundamentally 
different. 



1 1.4. Test of chromospheric radiation approximation 

There are no analytical tests that can be used to test our recipes 
of the chromospheric radiation described in Sect. 



8.2 The best 



we can do is to make a comparison with a simulation where the 
approximated processes have been calculated in detail. For this 
purpose we use ID radiation hydrodynamic simulations calcu- 
lated with the RADYN code, see e.g., Carlsson & Stein (1995, 
2002). This code solves the one-dimensional equations of mass, 
momentum, energy, and charge conservation together with the 
non-LTE radiative transfer and population rate equations, im- 
plicitly on an adaptive mesh. 

Since we have used such a simulation to determine the free 
parameters in our recipes, we use a simulation with a different 
velocity field for the test. Fig.|4]shows the comparison of the ra- 
diative cooling for one timestep in the RADYN simulation with 
the cooling calculated with the Bifrost recipes. The timestep 
shown has a strong wave that is close to steepening into a shock 
at lg(mc)=-4.2 (height of 1.2 Mm). The maximum temperature 
at the wave crest is 7000 K. Above the wave the temperature 
rises rapidly into the corona. At lg(mc)=-5.5 the temperature is 
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Fig. 3. The result of the Orszag-Tang test on a 256 x 256 grid, showing gas pressure (upper left), magnetic pressure (upper right), 
divergence of the velocity V • u (lower left) and the rotation of the velocity (V x u), (lower right). This figure can be directly 
compared with fig. 3 of Ryu et al. ( 1998| l. 



9000 K. The Bifrost approximations come close to describing 
the cooling in both hydrogen and calcium except that the cool- 
ing is overestimated in hydrogen just below the transition region 
(left part of the figure). Inspection of the cooling as function of 
time at a given height shows that the recipe for hydrogen typ- 
ically overestimates the cooling in the hot phases but does not 
include heating in the cool phases with the integral over time 
being close to the RADYN results. For further tests of the chro- 
mospheric radiation approximations see Carlsson & Leenaarts 
(in preparation). 



11.5. Full radiative transfer test in an isothermal scattering 
atmosphere 



The radiative transfer equation with coherent scattering has 
analytical solutions in the case of a static ID plane-parallel 
isothermal atmosphere if the photon destruction probability e 
is constant at all depths (e.g., Rybicki & Lightman 1979). The 
anisotropy of the radiation field needs to be restricted to linear 
dependence on the cosine of the zenith angle, fi - cos 6, mak- 
ing the Eddington approximation exact and reducing the prob- 
lem to solving a second-order equation for the mean intensity J. 
This setup is also known as the "two-stream" approximation. As 
second-order Gauss-Legendre quadrature yields an exact repre- 
sentation of integrals over a linear polynomial, a numerical result 
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Fig. 6. Numerical solution for the mean intensity 7nura as a function of optical depth t (left column), relative deviation A7rei from 
the analytical solution as a function of optical depth t (center column), and maximum relative deviation max(A7rei) as a function of 
iteration count (right column). The photon destruction probability e ranges from 10 ' (top row) to 10 (bottom row). The dashed 
and dot-dashed lines in the left column and center column mark the optical surface at t = 1 and the thermalization depth Ttherm- 
The dotted line in the center column indicates zero deviation. Dashed lines in the right column show the convergence speed for 
Gauss-Seidel corrections applied only during upsweeps, solid lines show the convergence speed for corrections appUed during both 
upsweeps and downsweeps. 



for J becomes directly comparable to an analytical solution (cf. 
Trujillo Bueno & Fabiani Bendicho|1995[ l. 



We set the radiative transfer solver to reproduce the mean 
radiation field 



An(T) = 



1 - 



1 + Ve 



(43) 
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Fig. 4. Radiative cooling as function of column mass for a de- 
tailed ID simulation with the RADYN code (solid) and with the 
Bifrost recipes for chromospheric radiation (dashed). Total cool- 
ing (black) and the contributions from hydrogen lines and the 
Lyman continuum (red) and lines from Call (blue). 




Fig. 5. Numerical solution for the mean intensity J in LTE (e - 
1.0) as a function of optical depth t. The dashed line at t = 1 
marks the optical surface. 



with the source function 



S(t) 



1 



B. 



(44) 



where t is the optical depth in the atmosphere, e is the photon 
destruction probability and B is the Planck function. The full ra- 
diative transfer module of the Bifrost code operates on a static 
3D atmosphere with a resolution of 50 x 50 x 120 grid points and 
with a horizontally homogeneous isothermal stratification; other 
physics modules are not used during the computation. The inter- 
polation algorithms needed for radiative transfer in 3D simula- 
tions are validated separately with the searchlight test described 
in |Hayek et al. (2010) . Optical depths are preset on a grid that 
is equidistant in logjQ t with 10 grid points per decade. Specific 
intensities are computed at /i = ±1/ V3 and arbitrary azimuth 
angles 0. The radiative transfer code uses double precision arith- 
metic to handle strong scattering at large optical depths, which 
appears due to the constancy of e in the atmosphere. The Planck 



function is set to Z? = 1 .0 in arbitrary units for all computations; 
it is also used as a first-guess source function for the solver 

In the LTE case (e = 1), the solver delivers I^ij) = 1.0 for 
outgoing intensities at /i = 1/ V3 and /"(t) = 1.0 - g^^^'" for 
ingoing intensities at fi - -1/ yf3 (see Fig.pll. The numerical so- 
lution is equivalent to the analytical solution given by Eq. 43 as 



the Gauss-Seidel solver uses the formal solution of the radiative 
transfer equation, which leads to identical expressions. In the 
scattering case, photon destruction probabilities assume values 
between e = 10"' (moderate scattering) and e = 10"^ (strong 
scattering), decreasing by factors of 10. The left column of Fig. 
|6]shows the numerical results for J. The mean intensity near the 
surface decreases for smaller e due to outward photon losses, 
and the thermalization depth Ttherm = 1 / ^ where the radiation 
field is completely thermalized (J ^ B), moves deeper into the 
atmosphere (dot-dashed line). At the smallest optical depths, the 
numerical solution delivers JnumiT - lO"**) ~ y/e for small e, as 



expected from Eq. 43 



The center column of Fig. |6] shows the optical depth- 
dependence of the total error of the numerical solution as the 
relative deviation between y,iuin and J-^^, 

A7rel(T) = — . (45) 

-/an(T) 

At large optical depths t » 1, radiative transfer is local and the 
total error depends on the source function gradient through the 
discretization error of the logarithmic optical depth grid. Using 
Eq. 44 one obtains the variation AS of the source function be- 
tween adjacent grid points with constant spacing A log r. 



AS 



A, 

T— — AlogT 
dr 



( 1 - Ve)B ( V3^Te" j A log r 

(1 - Vi)B(V3(T/Ttom)e-^^/^'^™)Alogi 



(46) 



In the LTE case, AS' = everywhere in the atmosphere since 
e = 1 .0 and S - B - 1 .0, independent of A log r, and the accu- 
racy of the numerical solution does not depend on the grid reso- 
lution. For e < 1 .0 and at optical depths t > Ttheim. the radiation 
field is entirely thermalized (ynum ~ B — \ .0), and A/iei vanishes 
since AS — > 0.0. The source function starts to decrease quickly 
through the decreasing mean intensity near Ttherm due to outward 
photon losses, causing a sharp increase in the error of the numer- 
ical solution. A/rei peaks in the translucent zone as AS is largest 
between 1.0 < t < T(heim- The magnitude of the peak grows 
with decreasing e through the (1 - ^|e) factor in Eq. 46 an upper 
limit is reached with e — > 0.0. It also follows from Eq. 46 that 
max(A7iei) scales with the grid resolution. Further up in the at- 
mosphere, the error decreases again through the finer grid spac- 
ing. At optical depths r <K Ttherm, the local discretization error 
becomes negligible due to the vanishing AS — > 0.0. However, 
as 7^ and I decouple from the local source function and the 
radiation field becomes anisotropic, the error becomes indepen- 
dent from the source function gradient. A/jei is dominated by the 
propagated error of outgoing radiation from deeper layers and 
is therefore constant. 

The right column of Fig. |6] shows the convergence speed of 
the numerical result to the analytical solution, measured through 
the maximum relative deviation max(Ayrei) in the atmosphere 
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Fig. 7. Ratio of the hydrogen level populations computed 
with the hydrogen ionization module to the populations from 
a detailed computation assuming statistical equilibrium. Upper 
panel: atomic hydrogen ground level; lower panel: proton den- 
sity. 



(see Eq. 45 i. Dashed lines represent computations where Gauss- 
Seidel corrections for the source function were applied only dur- 
ing upsweeps, while solid lines show the results when both up- 
sweeps and downsweeps were corrected, which increases con- 
vergence speeds by about a factor of two (see jTrujillo Buenoj 
|& Fabiani Bendicho|1995| l. Convergence is significantly slower 
when the thermalization depth moves to deep, very optically 
thick layers for small photon destruction probabilities, requiring 
about 250 iterations in the strongest scattering case. We find a 
numerical error of max(A7rei) ~ 3.7 • 10"^ at e = 10"*, similar to 



the error quoted in Trujillo Bueno & Fabiani Bendicho ( 1995 i. 



11.6. Hydrogen ionization test 

Non-equiUbrium hydrogen ionization generally produces atomic 
level populations that are neither in statistical equilibrium nor 
in Saha-Boltzmann equilibrium (LTE). However, in the case of 
simulations of the solar atmosphere spanning from the upper 
convection zone to the lower corona, two limiting cases are pro- 
duced: near the convection zone boundary the level populations 
obey LTE statistics; near the coronal boundary the populations 
obey statistical equilibrium because of the fast transition rates 
there. Both limiting cases are reproduced in statistical equilib- 
rium non-LTE radiative transfer codes. We will therefore com- 
pare the results from a Bifrost run with those of a statistical equi- 
librium code in those two limits. In the intermediate regime the 
statistical equilibrium is not valid and there will be large differ- 
ences between the two codes. 

We took a snapshot from a 2D simulation of the solar atmo- 
sphere that included the hydrogen ionization module (HION). 
This simulation had a grid size of 512 x 325 points, spanning 
from 1.5 Mm below the photosphere to 14 Mm above it. It in- 
cluded a weak magnetic field (average unsigned flux density in 
the photosphere of 0.3 G), thermal conduction and full, chromo- 
spheric and optically thin radiative transfer. The convection zone 
boundary was open, the coronal boundary used the methods of 
characteristics. 



We then computed the statistical equilibrium values of the 
hydrogen level populations from this snapshot using the radia- 
tive transfer code Multi3D treating each column as a ID plane- 
parallel atmosphere. As input we took the snapshot geometry, 
the mass density, the electron density and the temperature. We 
set the velocity field to zero. As model atom we used a 5- 
level plus continuum hydrogen atom with all radiative transitions 
from the ground level put in radiative equilibrium, as is assumed 
in the HION module. We treated all remaining lines assuming 
complete redistribution. 

Fig.jTJshows a comparison of the level populations obtained 
from the Bifrost and the Multi3D computation. Both the n = 1 
and the proton density are equal in the convection zone (below 
z = Mm), showing that Bifrost correctly reproduces LTE pop- 
ulations there. The proton densities in the corona are also equal 
for Bifrost and Multi3D. The n - \ populations in the corona 
are not identical. 

However, the differences are small, the populations are al- 
ways within a factor of two of each other The exact value of the 
coronal population density depends on the radiation field, and 
the simple HION assumption of a coronal radiation field that 
is constant in space does not capture the small-scale variations 
present in the Multi3D computation as indicated by the vertical 
stripes in the upper panel of Fig. [7] This striping is caused by 
variations in the atmospheric quantities down in the photosphere 
and chromosphere where the coronal radiation field is set. 

11.7. Boundary condition test 

The full 3D MHD equations allow a wide variety of wave modes 
that cross the boundaries at different angles depending on their 
origin and on the topology and strength of the magnetic field. A 
comprehensive test of all these modes falls outside the scope of 
this paper, but we will present a couple of tests to show the range 
of boundary condition behavior we have observed. The examples 
shown contain the full set of physics that Bifrost supports and 
thus are meant to represent the 'typical' production run the code 
is meant for, though with simpler geometries. Both of the tests 
shown were at the same time used to test the thermal conduction 
module as described in Sect. 111. 8] 

The first test is a ID model containing a photosphere, chro- 
mosphere, and corona that has a vertical extent z - 9 Mm, which 
is discretized on a grid containing 256 points with Az - 35 km 
throughout the computational domain. The model has solar grav- 
ity and a weak (0.1 G) vertical magnetic field. Radiative losses 
are included through the optically thin and chromospheric ap- 
proximations, but the full radiative transfer module is turned off. 
Thermal conduction is included and the corona is kept heated by 
maintaining the temperature at the upper boundary at 1 . 1 MK. 
The initial model is not perfectly in hydrostatic nor energetic 
equilibrium and the atmosphere responds by launching a set of 
acoustic disturbances at roughly the acoustic cut off frequency of 
3 minutes. These are initially of high amplitude, forming shock 
waves, but at later times damping out to much lower amplitudes 
and forming linear waves. In Fig.|8]we plot the vertical velocity 
M, as a function of height and time. Evidently, acoustic waves 
originate in the lower to middle chromosphere near z ~ 500 km 
and propagate upward steepening into shocks. At roughly 2 Mm 
they encounter the transition region before entering the corona 
proper. 

It is clear that while the boundary seems fairly well behaved, 
the strongest shocks do lead to some reflection, which we have 
measured to be of order 5% or less in terms of reflected energy. 
There may be several reasons for such reflection: One is that 
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Fig. 8. Vertical velocity m, in ID test model. The color scale is 
set to span +30 km/s with black representing upflow and white 
downflow. The chromosphere adjusts its structure, which ini- 
tially is not quite in hydrostatic equilibrium, by emitting acoustic 
waves at the cut-off frequency. These waves are initially strong 
enough to form shocks, as here during the first 1000 s of the 
simulation. 
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Fig. 9. Vertical velocity m- in 2D test model, red is downflow and 
blue upflow with a color scale set to ±50 km/s, as a function of 
time and height at the position x = 5 Mm in our 2D test model 
(see text). The disturbance clearly visible is a fast mode wave 
that originates due to a slight imbalance in the Lorentz force 
in the transition region in the initial state. The wave propagates 
fairly cleanly through the upper boundary, with a reflection of 
some 5%. 



high amplitude waves seem in general more difficult to transmit 
through a boundary. In addition, setting a hot plate as a tempera- 
ture boundary condition will give some reflection as the temper- 
ature is forced to a given value. 

The second experiment uses the same background atmo- 
sphere as in the ID test described above, but we have expanded 
to span 16.5 Mm in the horizontal direction forming a 2D at- 
mosphere that contains a magnetic field of greater complexity: 
This magnetic field is formed by inserting positive and negative 
magnetic polarities of ±1000 G strength at the bottom boundary, 
spanning 1 Mm and centered at x = 2.5 Mm and x - 13.5 Mm 
respectively, and then computing the potential field that arises 
from this distribution (see Fig. [TT] i. Note that the field is quite 
strong and that plasma beta is less than one in most of the mod- 
eled domain. Again there is a slight hydrostatic imbalance in 
the initial state and a transient wave is generated, in this case in 
the form of a fast mode wave originating close to the transition 
region. The Alfven speed is quite high, some 9 000 km/s near 
the transition region, falling to 1 000 km/s at the upper bound- 
ary, since much of the field has closed at lower heights in the 
atmosphere. In comparison, the speed of sound lies in the range 
100 km/s to 160 km/s in the corona. Thus, the generated fast 
mode, traveling essentially at the Alfven speed, propagates to the 
upper boundary very quickly, using only 3 s to cover the 7 Mm 
from the transition region to the upper boundary. The transient 
wave's amplitude on leaving the upper boundary depends on lo- 
cation but lies between -70 km/s and 150 km/s. In Fig. [9] we 
plot the vertical velocity as a function of height and time at the 
representative horizontal position x = 5 Mm. As in the previous 
example we do find a certain amount of reflection, but again find 
it to be energetically unimportant, at less than 5% of the energy 
contained in the original wave at all horizontal locations. Note 
that the reflected wave is re-reflected off the transition region 
and that this second wave seems very nicely transmitted through 
the upper boundary. 

Setting boundary conditions for the type of models discussed 
here is a compromise between long term stability and the best 



possible transmission of outwardly propagating modes. In the 
examples presented above we have attempted to show the results 
where stability is the most important factor, which is attained by 
limiting the velocity amplitude allowed in the ghost zones, as 
would be typical for a production run that aims to model the 
outer solar atmosphere for a duration of order several thousand 
seconds. This safety factor increases the amount of reflection 
seen, a compromise that allows us to run simulations for long 
periods of time. 

11.8. Thermal conduction test 

In order to test the thermal conduction module we present the 
temperature structure that arises in the two examples presented 
above in Sect. ll l.Tl In the first we consider a ID model with ver- 
tical magnetic field of the upper photosphere, chromosphere and 
corona in which the chromosphere is slowly relaxing by shed- 
ding acoustic waves at the cut-off frequency. The upper bound- 
ary temperature is set to 1 . 1 MK and conduction is the dominant 
term in the energy balance in the corona and transition region, 
down to temperatures of roughly 10 kK occurring at a height of 
z = 1 .2 Mm. After 3 hours solar time the amplitude of the acous- 
tic waves being generated is much reduced, and the location of 
the transition region moves by less than 300 km during a full 
wave period. In Fig. 10 we plot the temperature as a function of 



height for several timesteps during a wave period 3 hours solar 
time after the experiment began. Also plotted is the temperature 
profile that arises from a much higher resolution ID model com- 
puted with the TTRANZ code (Hansteen 1993) with the same 
upper boundary temperature. The latter model uses an adaptive 
grid (Dorfi & Drury 19871 that concentrates grid points in re- 
gions of strong gradients; in the present model the grid size is 
of order 80 m in the lower transition region. The temperature 
profiles in the two models are quite similar even though the grid 
spacing is much larger in the Bifrost run. Of course, not all as- 
pects of the relevant physics can be reproduced with Az - 35 km 



14 



Gudiksen et al.: The numerical code Bi frost 




Fig. 10. Comparison of high resolution ID model (solid line 
Az > 80 m) with Bifrost test run with Az - 35 km (dashed 
lines). The models are set up to have a temperature maximum of 
1 . 1 Mm at the upper boundary such that conduction dominates 
the energetics of the atmosphere in the corona and transition re- 
gion. The results from the Bifrost run are taken from the same 
model as used in Fig. |8] but at much later times (f > 8 000 s) 
when acoustic perturbations are largely damped and the ampli- 
tude of transition region motion is less than 300 km. 



in the high temperature gradient environment of the lower tran- 
sition region. We find sudden large temperature changes in tran- 
sition region grid points as the location of the transition region 
changes due to the passage of acoustic waves. These temperature 
changes become sources of (higher frequency) acoustic waves 
with amplitudes of order some 100 m/s that can be discerned on 
close inspection of Fig.|8] This artifact first disappears when the 
grid size is set to Az < 15 km (for typical coronal temperatures 
< 2 MK). 

In the second example we consider the case of rapid heating 
of coronal plasma in a magnetized atmosphere and follow how 
thermal conduction leads the deposited energy along the mag- 
netic field lines. We use the same atmosphere as described in the 
2D case given in the boundary condition test above. During the 
first second of the model run we deposit 50 J/m^ over a region 
spanning 100 x 100 km^ (an amount of energy equivalent to a 
large nano-flare) at location x - 1 Mm, z = 6.3 Mm, after which 
we allow the atmosphere to cool. At the upper boundary we set 
a zero temperature gradient boundary condition. 

Fig. 11 shows the temporal evolution of this model with 
s, and 10 s. The deposited energy 



snapshots taken at 0.6 s, 1. 
rapidly heats the coronal plasma, originally at 1 MK, to 2.4 MK 
at f = 1 s. The heat is efficiently conducted away from the site of 
energy deposition and has already after 0.6 s increased tempera- 
tures in a region several thousand kilometers along the magnetic 
field. After energy deposition ends the maximum temperature 
in the heated region decreases while the region itself continues 
to expand along the magnetic field towards the transition region 
and chromosphere. At f = 10 s the heated region has cooled 
quite a bit, but is still hotter than the ambient atmosphere and 
has spread out to form a loop-like structure. Note that the am- 
bient atmosphere is also cooling, the portions of the atmosphere 
connected via the magnetic field to the upper boundary cooling 
least rapidly. The width of the heated region is fairly constant. 




X [Mm] 

Fig. 11. Time evolution of 2D model which is heated for 1 s at 
position X = 7 Mm, z = 6.3 Mm with 50 i/nr' over a region span- 
ning lOOx 100 km^. The panels show the temperature at f = 0.6 s 
(top) when the maximum temperature is 2.25 MK, atf = 1.8 s 
(middle) when the maximum temperature is 1.63 MK, and at 
f = 10 s (bottom) where the maximum temperature is 1.16 MK. 
The temperature increases rapidly in the heated region, reach- 
ing 2.4 MK at 1 s, and plasma is heated by thermal conduction 
along the field as the plasma cools. Note that the upper boundary 
is set to have zero temperature gradient, so the entire atmosphere 
cools as well. Magnetic field lines are indicated with thin grey 
contours and contours of constant yS are shown with white num- 
bered lines. 



but some spreading perpendicular to the magnetic field is evi- 
dent in the last frame shown at f = 10 s. 



12. Parallelization 

Bifrost was written to be massively parallel. It employs the 
Message Passing Interface (MPl), because it is very well de- 
veloped and exists on almost all super computers. There are a 
number of other options for parallelization, including OpenMP 
and the thread () mechanism included in the different varia- 
tions of the C programming language, but we found MPl to suit 
our needs the best. 

MHD on a regular grid is trivial to parallelize by splitting the 
computational grid into subdomains and distribute one subdo- 
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Fig. 12. The relative number of computed gridpoints (solid) and 
the total data needed to be communicated (dashed) as a function 
of the number of internal grid points per dimension for a 3D 
cubic sub domain 



main to each node. In the case of Bifrost, computing the deriva- 
tives or interpolating the variables uses a stencil of 6 grid points. 
The two/three gridpoints nearest the edge of a subdomain, then 
relies on data that is outside the subdomain belonging to the lo- 
cal node, but instead belongs to the neighboring subdomain. The 
way this data is acquired most efficiently, depends on the prob- 
lem. The two possible solutions are to communicate the needed 
data to neighboring nodes every time it is needed, and the other 
is to supply each node with a number of "ghost cells" around 
its allocated subdomain so that the whole stencil used in for in- 
stance a derivative operator is present at the node that needs it. 
The extreme solution would be for each node to have a very large 
number of ghost cells, in the most extreme case the complete 
computational grid, and therefore never need to communicate, 
but that would make the code non-parallel. Choosing to include 
ghost cells around each sub domain makes it necessary to do 
more computations, since the ghost cells are copies of cells be- 
longing to neighboring nodes, but makes it possible to do less 
communication, as the ghost cells can be filled with the correct 
values from the neighboring nodes less often. The number of 
ghost cells chosen is therefore influenced by the relative impor- 
tance of the communication speed, the computation speed and 
the numerical scheme of the code. 

We have chosen to keep the communication between nodes 
to a minimum at the expense of doing more computations. Five 
ghost cells makes it possible to do both a derivative and an in- 
terpolation along the same coordinate direction without having 
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to communicate with neighboring nodes. That choice was made, 
because it allows us to get good performance even on machines 
with relatively slow internode communication speeds. It is very 
hard to quantify exactly when this choice is wise and when not, 
because so many parameters enter the problem. For instance, a 
global minimum operation scales with the number of nodes used 
and the communication speed, a simple neighbor communica- 
tion will depend both on the communication time, the physical 
setup of the nodes and the switches that connect them etc. But 
it is worth noting that typical communication times in large sys- 
tems are of millisecond scale, while the frequency of the cpus 
are in the gigahertz range. So from this very simple argument, 
it should be possible for a CPU to do roughly 10^ computations 
(multiplications, additions), in the time it takes to do one com- 
munication. 

Fig.[T2]shows how the data communicated increases with the 
number of internal gridpoints per node, and how the relative in- 
crease in computation decreases with the number of internal grid 
points. Both of the curves in Fig. [12] do not take into account 
communication speed or computation speed, so both curves can 
be shifted up and down the y-axis when applying them to a 
specific system. Communication can be split into initialization 
and actual communication of the data. In general the initializa- 
tion takes very long compared with the actual transmission of 
the data, so there is a large offset in the communication time, 
but the MHD part of Bifrost only communicates with neighbor- 
ing subdomains, so this offset depends on the dimensionality of 
the problem, and not on the number of cpus or internal points 
per subdomain. Since Bifrost was developed to handle as large 
computations as possible, the number of internal grid points will 
be rather high, so the rapid divergence between the two curves 
makes it plausible that doing the extra computations is the cor- 
rect choice. 

Scaling can be measured in two different ways, called strong 
and weak scaling. Strong scaling uses a set problem size and then 
timing measurements are made for different number of nodes, 
while weak scaling uses a set problem size per node, which 
means that when the number of nodes increases, the problem 
size also increases. Strong scaling is mainly a test of the com- 
munication overhead. If the communication takes up a constant 
amount of time, it should take a relatively larger and larger part 
of the run time as the number of cpus goes up. Weak scaling 
gives a measure of how well the code handles larger and larger 
problem sizes and number of cpus without being influenced by 
the raw communication time. 

Both strong and weak scaling tests of Bifrost have been 
performed on a number of computer architectures. We will 
here report on the results from a Cray XT4, with each node 
containing one quad-core AMD Opteron 2.3 GHz cpu and 
with a proprietary Cray Seastar2 interconnect, made available 
to us by the Partnership for Advanced Computing in Europe 
(PRACE), and from a Silicon Graphics ICE system, with each 
node containing two quad-core Intel Nehalem-EP 2.93 GHz 
cpus with Infiniband interconnect, located at NASA Advanced 
Supercomputing Division. 

Strong scaling tests were run with just the simplest con- 
figuration: pure MHD on an uneven staggered mesh, with an 
ideal EOS, without radiation, conduction or any other advanced 
physics or boundary conditions. The timing is performed on the 
Cray XT4 in such a way that the number of internal points per 
dimension is 30 or more, so the lower end of Fig. [12] is never 
reached. Such a test will show if there are communication bot- 
tlenecks in the code which will significantly slow the code down 
when running on a large number of nodes and when there is too 
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Fig. 13. Scaling when running a pure MHD case with 500^ grid- 
points on a Cray XT4 system with different number of cores (tri- 
angles), the theoretical scaling curve (dashed) and the theoretical 
scaling curve taking into account ghost cells (dotted) 



large a penalty due to the use of ghost cells. Fig. [13] shows the 
scaling results for a 500^ gridpoint run, and as the number of 
cores increases the relative amount of grid cells that are ghost 
cells increases, and consequently the code uses more and more 
time calculating the ghost cells compared to the internal cells. If 
the relative increase in ghost cells is taken into account, Bifrost 



scales very well. Fig. 13 also shows that if the number of cells 
per dimension for each core becomes less than 50, then effi- 
ciency of the code has dropped by about 35% compared to per- 
fect scaling, and consequently, each core should not get a com- 
putational sub domain which is smaller than 50 grid points on a 
side. 

The weak scaling tests were performed with a production 
setup of a solar simulation extending from the convection zone 
2.5 Mm below the photosphere to the corona 14.5 Mm above 
the photosphere. The horizontal extent was 24 x 24 Mm^. When 
the full radiation module was switched on in some of the weak 
scaling tests, scattering was included and the radiation field was 
described with 26 rays. The bottom boundary was transparent, 
the top boundary used the method of characteristics and hori- 
zontal boundaries were periodic. 

In theory Bifrost should be able to handle weak scaling very 
well, when only modules that use local data are used. If modules 
like full radiative transfer and Spitzer conductivity are included, 
the scaling is expected to drop below the theoretically best scal- 
ing, because these modules use non-local data and provide non- 
local results. When including modules using non-local data, the 
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Fig. 14. Weak scaling results for Bifrost on a Cray XT4 system 
when running MHD with a realistic EOS and chromospheric ra- 
diation (squares), MHD with a realistic EOS, chromospheric ra- 
diation and Spitzer heat conduction (diamonds) and MHD with 
a realistic EOS, chromospheric radiation, Spitzer heat conduc- 
tion and full radiative transport (triangles). Dashed lines show 
the average timing on the runs up to 256 cores. 
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Fig. 15. Weak scaling results for Bifrost on a Silicon Graphics 
ICE system when running MHD with a realistic EOS and chro- 
mospheric radiation (squares), MHD with a realistic EOS, chro- 
mospheric radiation and Spitzer heat conduction (diamonds) and 
MHD with a realistic EOS, chromospheric radiation, Spitzer 
heat conduction and full radiative transport (triangles). Dashed 
lines show the average timing on the runs up to 256 cores. 



drop in efficiency with larger number of cpus depends on the 
choice of, and implementation of, the solver as well as commu- 
nication time. For our tests the problem size was adjusted such 
that each core had a subdomain of 64 x 64 x 64 internal grid- 
points. The weak scaling results for Bifrost on the Cray XT4 are 
provided in Fig. 14 and show very little dependence on the num- 



ber of cores for each type of experiment. For the full radiative 
transfer case it is important to note that the number of iterations 
for convergence of the scattering radiative transfer problem de- 
pends on the change from the previous timestep, the aspect ratio 
between spacing in the vertical and horizontal directions and the 
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resolution. The large fluctuation in the timings for the full ra- 
diative transfer case seen in Fig. 14 is completely caused by the 
variations in the number of iterations per timestep for the dif- 
ferent cases and when this is taken into account the time per 
timestep does not increase with increased number of cores. Fig. 
15 shows the same experiments but run on the Silicon Graphics 
ICE system with Infiniband interconnect. There is again little de- 
pendence on the number of cores when the trend of increasing 
number of iterations with increased resolution for the radiation 
is taken into account. 



Figs. 14 15 also show that both the Spitzer heat conduction 



and full radiative transfer take up a large fraction of the time. 
Spitzer heat conduction increases the time by almost 50%, while 
full radiative transfer increases the computing time by more than 
a factor of 5. The large computational effort needed for the full 
radiative transfer is caused by the scattering iterations — typ- 
ically 3-15 iterations are needed for convergence within one 
timestep. In practical problems, the full radiative transfer does 
not have as large a penalty on the timing as it might seem, since 
the radiative transfer module only needs to be run in a fraction of 
the timesteps for chromospheric/coronal problems where scat- 
tering is important. In a typical solar simulation with 32 km 
horizontal resolution (the 1728 core points in Figs. 14 15 1 the 



timestep set by the Courant-Friedrich-Levy condition and the 
Alfven speed in the corona is 3 ms while the radiation is up- 
dated only every 300 ms which means that the radiative trans- 
fer increases the computing time only by 2%. For photospheric 
problems, the radiative timescale is comparable with the hydro- 
dynamic one and the radiative transfer module needs to be called 
for every timestep. On the other hand, there is no need to iterate 
when scattering is unimportant. For such simulations the radia- 
tive transfer increases the computing time by about 60% com- 
pared with the pure MHD case. 

13. Conclusion 



The development of numerical methods and computer power has 
made numerical simulations of stellar atmospheres highly rele- 
vant. It is now possible to create observational predictions from 
advanced realistic numerical simulations, and observations can 
therefore partially validate the results of the numerical simula- 
tions. If such predictions are made and confirmed by observa- 
tions it is likely that other predictions from the numerical sim- 
ulations are also correct, making it possible to get much more 
information about the stellar atmospheres than would be possi- 
ble through observations alone. To provide as good simulation 
results as possible, it is necessary to have a highly efficient and 
parallel numerical code. The tendency for modern super com- 
puters to be distributed memory systems, makes it extremely im- 
portant that numerical codes are highly parallel. For pure MHD 
that is not difficult to attain because MHD is a local process, but 
it becomes much more complicated for non-local processes or 
numerical solvers. 

The numerical MHD code Bifrost has been created to meet 
the requirements of present supercomputers. It is developed by 
a group of researchers, post docs and PhD students and provides 
a simple interface to include further developments in boundary 
conditions and physical regimes that are not included in the sim- 
plified core of the code. Several extension modules are already 
provided and several more are under development. Results using 
Bifrost have already been publis hed ([Hayek et al. 2010[|Carlsson| 
et al. 2010 ; Hansteen et al.|2010[|Leenaarts et al. In press, 201 1 



The pure MHD module of Bifrost has been extensively tested 
through standard tests and has performed well. There are at the 



moment eight individual modules that have been finished and 
tested, and several more are under development and testing. The 
finished modules include a number of equations of state, radi- 
ation transport and thermal conduction and tests of them have 
been presented and all produce results according to expectations. 

Bifrost has been tested for scaling performance. Both strong 
scaling and weak scaling results are very good on the two sys- 
tems we have tested on. We would ideally have liked to test on 
a system where the interconnect is slower than on the Silicon 
Graphics ICE and Cray XT4 systems to get further knowledge 
about the bottlenecks a slow interconnect would present for 
Bifrost. The hardware communication architecture plays a role 
on the scaling behavior of Bifrost, but these would most likely 
be more severe if we had made a choice of doing more commu- 
nication. Since Bifrost is primarily designed to run large simu- 
lations, using a large number of cores, we believe we have made 
the correct prioritization in choosing more computations over 
communication. 

The very good scaUng performance and the modules already 
developed will make it possible to simulate the whole solar at- 
mosphere from the top of the convection zone to the corona with 
a degree of realism that has not been attained before. It has be- 
come more and more clear that the solar atmosphere cannot be 
split into the traditional separate layers, the photosphere, chro- 
mosphere, transition region and corona. The solar atmosphere is 
one large connected system, and it is necessary to include the 
whole atmosphere to attain credible results. The consequence is 
that the code used for such a simulation will have to include 
the special physics important in each layer, making the numeri- 
cal code much more complex than a numerical code designed to 
deal with just one of the layers. Bifrost is uniquely qualified for 
that task. 

The relative ease of creating new setups for simulations, in- 
clusion of special physics and boundary conditions, make it pos- 
sible to use Bifrost for detailed solar atmosphere simulations and 
for stellar atmospheres in general. Several investigations using 
Bifrost have been done or are under way, several of these in- 
cluding PhD students who have been able to use this 'state of the 
art' numerical code with relatively little instruction. There are a 
number of modules being developed, including a module that 
can follow the ionization states of heavy elements and a module 
that introduces a modified Ohm's law. There is a trend towards 
computer systems using the large raw floating point performance 
of Grapical Processing Units (GPUs), and a new parallelization 
module for such systems is also under development. 

The very good scaling performance of Bifrost makes it pos- 
sible to make simulations of stellar atmospheres with a very 
large resolution while still encompassing a large enough volume 
to make the simulation realistic, and include a large number of 
physical effects. The results can be used to predict observational 
effects, which might earlier have lead to wrong diagnostics of 
the physical parameters in the solar atmosphere. 



Acknowledgements. We would like to thank the referee for very useful com- 
ments, that helped us make improvements. We would like to thank Robert F. 
Stein, Andrew McMurry and Colin Rosenthal for working out the formalism 
for the characteristic boundary conditions. This research was supported by the 
Research Council of Norway through the grant "Solar Atmospheric Modelling" 
and through grants of coinputing time from the Programine for Supercomputing 
and through grants SMD-07-0434, SMD-08-0743, SMD-09-1128, SMD-09- 
1336, and SMD-10-1622 from the High End Computing (HEC) division of 
NASA. The authors wish to thank PRACE for opportunity to run experiments on 
HPC centers in Europe. J.L. acknowledges financial support from the European 
Commission through the SOLAIRE Network (MTRN-CT-2006-035484) and 
from the Netherlands Organization for Scientific Research (NWO). 



18 



Gudiksen et al.: The numerical code Bifrost 



References 

Abbett, W. P. 2007, ApJ, 665, 1469 

Arber, T. D., Haynes, M., & Leake, J. E. 2007, ApJ, 666, 541 
Arnaud, M. & Rothenflug, R. 1985, A&AS, 60, 425 

Carlsson, M., Hansteen, V. H., & Gudiksen, B. V. 2010, 

Mem. Soc. Astron. Italiana, 81, 582 
Carlsson, M. & Stein, R. F. 1995, ApJ, 440, L29 
Carlsson, M. & Stein, R. F. 2002, ApJ, 572, 626 

Carlsson, M., Stein, R. F, Nordlund, A., & Scharmer, G. B. 2004, ApJ, 610, 
L137 

Cheung, M. C. M., Schiissler, M., & Moreno-Insertis, F. 2007, A&A, 467, 703 
De Pontieu, B., Hansteen, V. H., Rouppe van der Voort, L., van Noort, M., & 

Carlsson, M. 2007, ApJ, 655, 624 
DeVore, C. R. 1991, Journal of Computational Physics, 92, 142 
Dorfi, E. A. & Drury, L. O. 1987, Journal of Computational Physics, 69, 175 
Galsgaard, K. & Nordlund, A. 1996, J. Geophys. Res., 101, 13445 
Gudiksen, B. V. & Nordlund, A. 2005, ApJ, 618, 1020 
Gustafsson, B., Bell, R. A., Eriksson, K., & Nordlund, A. 1975, A&A, 42, 407 
Hansteen, V. 1993, ApJ, 402, 741 

Hansteen, V. H. 2004, in lAU Symposium, Vol. 223, Multi-Wavelength 
Investigations of Solar Activity, ed. A. V. Stepanov, E. E. Benevolenskaya, 
& A. G. Kosovichev (Cambridge University Press: Cambridge), 385 
Hansteen, V. H., Carlsson, M., & Gudiksen, B. 2007, in Astronomical Society 
of the Pacific Conference Series, Vol. 368, The Physics of Chromospheric 
Plasmas, ed. P. Heinzel, I. DorotoviC, & R. J. Rutten, 107 
Hansteen, V. H., De Pontieu, B., Rouppe van der Voort, L., van Noort, M., & 

Carlsson, M. 2006, ApJ, 647, L73 
Hansteen, V. H., Hara, H., De Pontieu, B., & Carlsson, M. 2010, ApJ, 718, 1070 
Hayek, W., Asplund, M., Carlsson, M., et al. 2010, A&A, 517, A49 
Heggland, L., De Pontieu, B., & Hansteen, V. H. 2007, ApJ, 666, 1277 
Heinemann, T, Nordlund, A., Scharmer, G. B., & Spruit, H. C. 2007, ApJ, 669, 
1390 

Hyman, J. M. 1979, in Proceedings of the third IMACS international symposium 
on computer methods for partial difierential equations, ed. R. Vichnevetsky, 
R. Stepleman, Vol. 3, International Association for Mathematics and 
Computers in Simulation (Dept. of Computer Science, Rutgers University, 
New Brunswick, N.J 08903 USA: IMACS), 313 
Isobe, H., Proctor, M. R. E., & Weiss, N. O. 2008, ApJ, 679, L57 
Judge, R G. & Meisner, R. W. 1994, in ESA Special PubUcation, Vol. 373, 
Solar Dynamic Phenomena and Solar Wind Consequences, the Third SOHO 
Workshop, ed. J. J. Hunt, 67 
Keller, C. U., Schussler, M., Vogler, A., & Zakharov, V. 2004, ApJ, 607, L59 
Leenaarts, J. & Carlsson, M. 2009, in Astronomical Society of the Pacific 
Conference Series, Vol. 415, Astronomical Society of the Pacific Conference 
Series, ed. B. Lites, M. Cheung, T. Magara, J. Mariska, & K. Reeves, 87 
Leenaarts, J., Carlsson, M., Hansteen, V., & Gudiksen, B. In press, 2011, A&A 
Leenaarts, J., Carlsson, M., Hansteen, V., & Rutten, R. J. 2007, A&A, 473, 625 
Martfnez-Sykora, J., Hansteen, V, & Carlsson, M. 2008, ApJ, 679, 871 
Martmez-Sykora, J., Hansteen, V, & Carlsson, M. 2009a, ApJ, 702, 129 
Martfnez-Sykora, J., Hansteen, V, DePontieu, B., & Carlsson, M. 2009b, ApJ, 
701, 1569 

Martrnez-Sykora, J., Hansteen, V., & Moreno-Insertis, F. 2010, ArXiv e-prints 
Miyama, S. M., Tomisaka, K., & Hanawa, T. 1999, Astrophysics and Space 

Science Library, Vol. 240, Numerical Astrophysics (Springer, Dordrecht) 
Nordlund, A. 1982, A&A, 107, 1 

Nordlund, A. & Galsgaard, K. 1995, A 3D MHD Code for Parallel Computers, 
Tech. rep.. Astronomical Observatory, Copenhagen University 

Orszag, S. A. & Tang, C.-M. 1979, Joimial of Fluid Mechanics Digital Archive, 
90, 129 

Press, W. H., Teukolsky, S. A., Vetterling, W. T, & Flannery, B. P 1992, 
Numerical recipes in FORTRAN. The art of scientific computing (Cambridge 
University Press, Cambridge) 

Rempel, M., Schiissler, M., Cameron, R. H., & Knolker, M. 2009, Science, 325, 
171 

Rybicki, G. B. & Lightman, A. P. 1979, Radiative processes in astrophysics 

(Wiley-Interscience, New York) 
Ryu, D., Mmiati, P., Jones, T. W, & Frank, A. 1998, ApJ, 509, 244 
Schafi'enberger, W, Wedemeyer-Bohm, S., Steiner, O., & Freytag, B. 2005, in 

ESA Special Publication, Vol. 596, Chromospheric and Coronal Magnetic 

Fields, ed. D. E. Innes, A. Lagg, & S. A. Solanki 
Shull, J. M. & van Steenberg, M. 1982, ApJS, 48, 95 
Skartlien, R. 2000, ApJ, 536, 465 
Sod, G. A. 1978, Journal of Computational Physics, 27, 1 
Solium, E. 1999, Master's thesis. University of Oslo 

Spitzer, L. 1956, Physics of FuUy Ionized Gases (New York: Interscience 
Publishers, New York) 



Stein, R. F., Lagerfjard, A., Nordlund, A., & Georgobiani, D. 2010, in American 
Astronomical Society Meeting Abstracts, Vol. 216, American Astronomical 
Society Meeting Abstracts, 2 1 1 .03 

Stein, R. E & Nordlund, A. 2006, ApJ, 642, 1246 

Stein, R. E, Nordlund, A., Georgoviani, D., Benson, D., & Schaffenberger, 
W. 2009, in Astronomical Society of the Pacific Conference Series, Vol. 
416, Astronomical Society of the Pacific Conference Series, ed. M. Dikpati, 
T. Arentoft, 1. Gonzalez Hernandez, C. Lindsey, & E. Hill, 421 
Thompson, K. W. 1987, Journal of Computational Physics, 68, 1 
Tortosa-Andreu, A. & MorenoJnsertis, E 2009, A&A, 507, 949 
Toth, T. & Odstrcil, D. 1996, Journal of Computational Physics, 128, 82 
Trujillo Bueno, J. & Eabiani Bendicho, R 1995, ApJ, 455, 646 
Tsuji,T. 1973, A&A, 23,411 

Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635 
Vogler, A., Shelyag, S., Schussler, M., et al. 2005, A&A, 429, 335 
Wedemeyer-Bohm, S. & Steffen, M. 2007, A&A, 462, L31 
Woodall, J., Agundez, M., Markwick-Kemper, A. J., & MUlar, T. J. 2007, A&A, 
466, 1197 



Appendix A: Characteristic boundary conditions 

The boundary equations in terms of the primitive variables 
ip, u, e, B) can be written in the following form 



dp 

dt 



~dt 



- 1^1 d2 + ^(ds + de) + ^(dj + Jg) 
-(uh • Vh)p - pVh • uh (A. 1) 



Ryi-ds + d4) + ^-^R^i-ds + 4) 

ct 



dUy 

~dt 



iHp+w) 1 

. -(uh-VhK+ (Bh-Vh)5;, (A.2) 

p ox jUoP 

Rxidi - d4) + ^^Ryi-ds + de) 



+ ^-^Ryidj-ds) 



iHp+w) 1 

. -(uh-VhK + (BH-VH)fiy (A.3) 



dt 



2c: 



de 
dt 



dB, 
dt 



cl 



[c+a+{di - de) + c-a-{di + d^)] 



-(Uh • Vh)m^ +g+ (Bh • Vh)Bz 

y^op 



-(|-| +\{e + P)a+{d5 + a^d6) 



+ ]^{e + P)a^{di + d^) 



-Vh ■ (eufl) - /^Vh • Uh + e 



Ry{d^ + di) + —Rxids + de) 



— -RJd7 + da) 
-(uh • Vh)B;, + B^iWn ■ Uh) + (Bh • Vh)mx 



(A.4) 
(A.5) 



(A.6) 



(A.7) 



dBy 



dt 



-RAd3 + cIa) + — Ry{ds + df,) 
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(a- dP 



a+ 



-(uh • VH)By + ByiVu ■ Uh) + (Bh • V^)uy 
-(uh • Vh)Bz - Bz(Vh • Uh) + (Bh • Vn)u, 



where we have defined the quantities 

S7 



- sign(fiz) 



Rr = 



al = 



Bh 

2 2 
c: - ct 



R. = 



Bh 



2 2 

- ~ 2 2 

C+ - ct 



(A.8) 
(A.9) 



d-i — (Mj + C-) 
du 



+C-a- 

d% = (mj-c-) 
du 



p dz 

RxCsa+ dBy, 

dz ^/^oP dz 
dP 



+ s^R^c^a^^^ + s^RyC+a+ 



du, 
' dz 
RyCsa+ dBy 
-y^uop dz 

dUr 



dUy 

di 



duv 



-C_Qr_ 



- s^R_,c^a+ — s.RyC+a+ - 

p dz dz ' dz 

RxCf.a+ dB, RyCsa^ dBy 



dz 



JiM)P dz 



\/PqP dz 
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(A.18) 



(A. 19) 



Along inflowing characteristics the characteristic derivatives 
d are changed to provide transmitting boundaries by requiring 
that the incoming characteristics remaining constant. Thus, we 
require that the boundary conditions for the incoming character- 
istics must satisfy the static case where u = and ^ = 0, so 
that 



using the velocities 



d = 



dP 
dp 

J"0P 



e + P IdP 



cl = 



cl + cl 



cl + cl 



- etc, 



2^2 



Note the unorthodox ordering of the characteristics: by conven- 
tion these are usually ordered by the amplitude of the character- 
istic speeds X-,, this was not known by us at the time these equa- 
tions were derived. Therefore, the characteristic z derivatives d 
are given by 



d\ — 



dB^ 



d2 
d3 



Uzie + P) 



dp de 



dz " 

{Uz + cM-s.Ry 



dz 

du,, 
Ik 



dUy 



Ry dBj, 



Rr dB, 



dz y/fioP dz 



(u^-Cz)\szRy 



dux 

dz 



SzRx 



dUy 

Ik 



(A.IO) 
(A.11) 
(A.12) 

(A.13) 

(A. 14) 



+ - 



Ry dB,, R^ dBy 



(Mz + c+) 
du 



TioP dz 

a. dP 



Jmp dz 



s^RjtC^a^—— 
p dz dz 



dUy 
^ dz 



+c+a+ 



dz 



RxCsO- dBx ^ RyC^a^ dBy 



Jpop dz 



JpoP dz 



(A.16) 



a+ dP 



dUy 



dUy 



de = (mz - c+) — -T- -I- s^Rj,c-a- — + sJtyC-U- — 
^ p dz dz dz 



-c+a+ 



duz 
d^ 



RxCsO- dBx ^ RyC^a- dBy 



JjloP dz 



JjloP dz 



di = (S-iC("=°'),- . 



(A.20) 



The components of C*"""^ are for the magnetic field equation 
and the density equation. For the other equations: 



r(«=o) 

(-(u=0) 
^(u=0) 



Q 



1 _a 

p dx 

1 d 



P + 



pdy 



B^ 



+ (Bh • Vt,)B, 

f^oP 

= --^1^+:^— 1 + — (Bh-Vh)Bv 



g + (Bh • Vh)fi, . 

f^OP 



Assuming 2 = 0, the non-zero incoming characteristic deriva- 
tive vector d can be calculated to be: 



4 = - 



Bi 



-(Rh X Vh) P + + — (Vh X Bh) 

2jUo/ j"o 



4 = 



c^a+\g+ (Bh- Vh)Bz 

Pop 



+c_a_^(Rh-Vh)|P + 
P 



dL = -4 



4 = 



4 



c-a-\g+ (Bh- Vh)fiz 



-c^a^^(Rh-Vh)|P+;^) 



-4 



(A.21) 

(A.22) 
(A.23) 

(A.24) 
(A.25) 
(A.26) 

(A.27) 

(A.28) 



. Then either di or d'. is chosen depending on the sign of Ai at 
the boundary; i.e. whether or not the characteristic is in- or out- 
flowing. 



Appendix B: Time-dependent H2 formation 

The module that computes non-equihbrium hydrogen ioniza- 
tion has been extended to include a ninth equation for time- 
dependent H2 formation. It is given by: 



(A.17) F9 = ^ - ^ (CsHWi - Cmmmni) -1=0, 



Af 



(B.l) 
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with «°, the H2 population of the previous timestep, At the 
timestep and «i the population of atomic hydrogen in the ground 
state. The rate coefficients are given by 

C3H = . (B.3) 

The values for a, /3 and y are taken from the UMIST database 
( |Woodall et al.||2007| iwww.udf a.net ); K(T) is the chemical 
equilibrium constant taken from |Tsuji| ( |1973| l. The derivatives 
of the functional Fg with respect to the dependent variables are 

dFg At I dCiii 3 5Ch2h \ ,„ 

-"1 T^^nmni , (B.4) 



dT n°2 \ 5r ' dT 

dFg _ At 
drii 

dFg _ 1 At 
dnm n° n° 



— (SCsH"! - Cmunm) , (B.5) 
Cmnni . (B.6) 



The rate equation for rii from Leenaarts et al. ( 2007| l is mod 



ified to include source and sink terms due to H2: 

f 6 



ni At 
F3 = — 

n° n° 



^ rijPji + Cmnnmni - C^nnl 



1=0, (B.7) 



with n° the ground state hydrogen population from the previous 
timestep. The derivatives of this equation and equations express- 
ing energy conservation and hydrogen nucleus conservation (see 
Leenaarts et al.|2007) l are modified correspondingly. 



